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Abstract. New axisymmetric (2D) models for dust-driven winds of C-stars are presented which include hydrody- 
namics with radiation pressure on dust, equilibrium chemistry and time-dependent dust formation with coupled 
grey Monte Carlo radiative transfer. Considering the most simple case without stellar pulsation (hydrostatic inner 
boundary condition) these models reveal a more complex picture of the dust formation and wind acceleration 
as compared to earlier published spherically symmetric (ID) models. The so-called exterior k- mechanism causes 
radial oscillations with short phases of active dust formation between longer phases without appreciable dust 
formation, just like in the ID models. However, in 2D geometry, the oscillations can be out-of-phase at different 
places above the stellar atmosphere which result in the formation of dust arcs or smaller caps that only occupy 
a certain fraction of the total solid angle. These dust structures are accelerated outward by radiation pressure, 
expanding radially and tangentially like mushroom clouds, while dust-poor matter is falling back towards the star 
at other places. A highly dynamical and turbulent dust formation zone is created in this way, which again leads to 
inhomogeneous dust production. Further away from the star, flow instabilities (e. g. Rayleigh- Taylor) have time 
to fragment the outward moving arcs and shells to produce numerous small-scale cloud-like sub-structures. 
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1. Introduction 



Asymptotic Giant Branch (AGB) stars are known to gen- 
erate massive stellar winds with mass loss rates up to 
lCT 4 M /yr (Habing 1996, LeBertre 1997, Wallerstein & 
Knapp 1998), which drives them towards the planetary 
nebula phase. According to the currently available spheri- 
cally symmetric (ID) model calculations (Winters et al. 
2000, Homer et al. 2003, Sandin & Homer 2003, Schirr- 
macher et al. 2003), these winds are driven by a combi- 
nation of stellar pulsation and radiation pressure on dust 
grains. An exterior K-mechanism leads to the production 
of radial dust shells in more or less regular long time in- 
tervals, even if the pulsation of the star is neglected in the 
model (Fleischer et al. 1995, Hofner et al. 1995). 

The formation of radial dust shells, however, does 
not seem to be in general agreement with recent in- 
terferometric observations. Infrared speckle observations 
have clearly indicated that the innermost dust distribu- 
tion around red giants can deviate strongly from spher- 
ical symmetry (Weigelt et al. 2002, Monnier et al. 2004, 
Tuthill et al. 2005) with evidence for moving clouds, arcs 
and "bars". The new IR Vlti instruments Amber and 
Midi will allow for an even more detailed view on the 
dust formation and wind acceleration zones of AGB stars 
with the ability to detect wind asymmetries and cloud- 
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like inhomogeneities, and to track individual dust shells 
or clouds. 

These observations call for likewise detailed 3D mod- 
els. Without such models, the interpretation of the new 
interferometric data will necessarily remain on a more 
or less phenomenological level which is unsatisfactory be- 
cause it leaves most of the important physical questions 
unanswered. 

The major aim of this work is to advance from ID 
to 2D dust-driven wind models. These models will enable 
us to discuss whether dust-driven winds remain spheri- 
cally symmetric (if the initial and boundary conditions 
are spherically symmetric) or whether instabilities lead 
to spontaneous symmetry breaking, structure formation 
and more complex flow patterns, which might explain the 
observed irregular dust structures. Post-processing of the 
model results will allow us to calculate simulated images 
and visibilities (see Woitke & Quirrenbach 2005) which 
can be directly compared to observations. 



2. The model 

2.1. The equations 

The dust-driven stellar wind is described by the equation 
of continuity (Eq.QJ, the equation of motion (Eq.[2J), the 
gas energy equation (Eq.|2J), a system of moment equa- 



tions describing the timc-dcpcndcnt formation of dust par- 
ticles from the gas phase (Eq.|3J , and the radiative transfer 
equation (Eq.EJ: 
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(pv) + S7(pvo V + p) = p(a rad - ag rav ) (2) 

pe) + V([pe +p]v) = pV-(a r ad-a g rav) + P<3rad (3) 



^(pLi) + V( P L jV ) = Vi j / 3 J± + I X nct Ph-x (4) 
n-VI v {r,n) = ry„(r,n) - < xt (r) I v {r,n) (5) 

p is the mass density, e = ei n t + eki n the internal (including 
kinetic) gas energy and p the thermal gas pressure. a grav 
is the gravitational deceleration and a ra( j the acceleration 
by radiation pressure on dust and gas. Q r ad is the net 
energy exchange rate of the gas with the radiation field 
(negative for radiative cooling). For further explanations 
of the introduced quantities see Sects . 1231 and l2~4l 

The dust/gas mixture is approximated by a sin- 
gle hydrodynamic fluid with bulk velocity v (neglect- 
ing drift velocities 1 ). However, since dust and gas are 
known to decouple thermally from each other easily un- 
der the predominant conditions in red giant winds (e. g. 
Woitke et al. 1999), two distinct temperatures are intro- 
duced, namely the gas temperature T g and the dust tem- 
perature Td. Whereas T g is a result of the solution of the 
gas energy equation (Eq.|2|) , we assume that the dust tem- 
perature Td relaxes instantaneously to its local radiative 
equilibrium value, which is part of the results of our ra- 
diative transfer treatment (see Sect. 12. 4|) . 



2.2. Hydrodynamics 

The multi-dimensional models are developed in the 
frame of the Flash 2.4 hydrocode (Fryxell et al. 2000), 
which is an explicit, finite volume, high-order Godunov- 
type hydro-solver which uses Adaptive Mesh Refinement 
(AMR) and is fully parallel (MPI). Spherical coordinates 
(r, 9, (f>) are used either in ID (spherically symmetric mod- 
els) or 2D (axisymmctric models). 

The source terms appearing on the r.h.s. of the hydro- 
equations (Eqs. 1 to 4) are treated by operator splitting, 
where the time-integration is carried out by means of the 
implicit extrapolation solver Limex4.2A1 (Deuflhard & 
Nowak 1987). 

Gravity and radiation pressure are calculated by 



Ctrad 



t£* GM r 



47rr 2 c 



(6) 



where e r is the unit vector in the radial direction, G 
the constant of gravity, the stellar luminosity, K e xt 
the grey (dust + gas) extinction coefficient per mass and 

1 See Simis et al. (2001) and Sandin & Hofner (2003) for the 
effects of dust drift velocities in ID models. 



c the speed ot light. Ihe enclosed mass is approximated 
by M r ^iM r (t — Q), i.e. we assume a constant spherically 
symmetric gravitational potential as calculated for the ini- 
tial state. The radiation pressure is assumed to be strictly 
directed into the radial direction. 

Net radiative heating of the gas is calculated by 



Q r ad-4^f b a s s (T r 4 ad -T g 4 



(7) 



assuming LTE, where Jcf^s is the gas absorption coefficient 
per mass and T ra d is a measure for the mean intensity 
(J=^T r \ d , seeSect.Q 

As equation of state, we apply the ideal gas law for 
atomic hydrogen with 7 = 5/3 (disregarding H2), i.e. 
p = 4 pT g and e mt — ^zjp/p, where k is the Boltzmann 
constant and p — 1.28 amu the mean molecular mass. 

2.3. Chemistry and dust formation 

The time-dependent formation of dust in the carbon-rich 
case is calculated according to the moment method de- 
veloped by Gail & Sedlmayr (1988), Gauger et al. (1990) 
and Dominik et al. (1993). Assuming spherical dust par- 
ticles, the moments of the dust size distribution function 
f(V,r,t) are defined as 

/>oo 

pL j (r,t)= f(V,r,t)V^ 3 dV (j e {0, 1, 2, 3}) , (8) 



where V is the dust particle volume and Vi the minimum 
volume of a large molecule to be counted as dust parti- 
cle. The temporal change of these dust moments, as af- 
fected by nucleation, growth and evaporation, are calcu- 
lated according to Eq. (4). Here, J* [cm~ 3 s^ 1 ] is the nucle- 
ation rate (i. e. the creation rate of seed particles) which 
is calculated for pure carbon chain molecules according to 
Gail et al. (1984). x nct [cm/s] is the net growth velocity of 
already existing dust surfaces by the accretion and ther- 
mal evaporation of carbon and hydrocarbon molecules. 
The functional dependencies of these quantities are 
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The various particle densities n k entering into this de- 
scription are calculated by assuming chemical equilibrium 
in the gas phase among H, H 2 , C, O, CO, H 2 0, C2, C3, 
CH, C2H, C2H2, CH4 according to the local mass density 
p, the local gas temperature T g , and the local carbon ele- 
ment abundance ec in the gas phase. The latter obeys the 
following conservation law 



ec + 



1.427 amu 



Vn 



L3 = e c = const 



(12) 



Vo = 8.78 • 10 -24 cm 3 is the monomer volume of graphite 
and 1.427 amu the conversion factor between mass density 
p and hydrogen nuclei density nm\ for solar abundances. 
All other element abundances than carbon are assumed 



to have solar values. We take the oxygen abundance eo = 
10 8.87-i2 f rom (Grevesse & Noels 1993) and consider the 
carbon-to-oxygen ratio C/0 = ec/eo as a free parameter. 

Dust evaporation is only included in terms of negative 
X nct in case of under-saturation. Negative particle fluxes 
in size space through the lower integration boundary at 
Vi are neglected. For more details, see Gail & Sedlmayr 
(1988), Gauger et al. (1990) and Dominik et al. (1993) 

2.4. Radiative transfer 

Radiative transfer in dusty winds is a complicated prob- 
lem. The optical depths are typically of the order of unity 
and the opaque regions can be confined into shells or 
clouds which move in an otherwise transparent medium. 
The opaque dust configurations scatter, re-emit and cast 
shadows, thereby illuminating the optically thin regions in 
between in complicated ways. Finding a efficient radiative 
transfer method to resolve these problems is the key to 
develop multi-dimensional models of dust-driven winds. 

The absorption and thermal re-emission of radiation 
by dust grains leads to a very fast relaxation of the in- 
ternal dust temperature to the local radiation field with 
cooling timescales of the order of milliseconds for small 
amorphous carbon grains (Woitke et al. 1999). This fast 
relaxation introduces a stiff coupling between hydrody- 
namics and radiative transfer. Any explicit scheme that is 
based on formal solutions of the radiative transfer problem 
(e. g. short characteristics method) has the disadvantage 
to slow down the computational timestep (typically a few 
10 4 seconds) to this radiative cooling timescale, which is 
not acceptable. 

We therefore seek for a numerical method that is ca- 
pable to solve the continuum radiative transfer problem 
under the auxiliary condition that the dust component is 
in radiative equilibrium 

J 4lfj„ du = J 4lfB v {T d ) dv , (13) 

which de-stiffens the physical problem by the elimina- 
tion of the shortest characteristic timescale via assuming 
a quasi-equilibrium. K^bs* i s the dust absorption coeffi- 
cient, J v (r) — f I v (r,n) d 2 n the mean spectral inten- 
sity, I v {r, n) the spectral intensity in direction of the unit 
vector n and B v the Planck function. In 2D or 3D, state- 
of-the-art Monte Carlo methods can do this job with the 
same accuracy as ray-based methods at a similar level of 
computational costs (e.g. Pascucci et al. 2004). 

We have therefore developed a new software module 
for the FLASH-code which solves the LTE continuum ra- 
diative transfer with radiative equilibrium in spherical co- 
ordinates by a Monte Carlo (MC) method. According to 
our knowledge, this is actually the first report of such an 
approach in multi-dimensional fluid dynamics. 

We have implemented the MC method of Niccolini, 
Woitke & Lopez (2003) which solves the frequency- 
dependent radiative transfer problem with angle-depen- 
dent scattering. However, for this paper, we will only con- 



sider the grey case with isotropic scattering and simple 
frequency-averaged opacities, because we want to validate 
our new implementation by comparing the results of ID 
spherically symmetric winds to former publications, where 
these opacities have been used (Fleischer et al. 1995, 
Hofner et al. 1995): 

cxt I •£ cxt / 1 A \ 

^cxt — ftg as "I" /-Must 

= 2.10- 4 cm 2 g- 1 (15) 
«dSt = \Q\T d )L 3 with Q'(T d ) = 5.9 • T d (16) 

The constant gas opacity is taken from Bowen (1988) and 
the Rosseland mean opacity for small amorphous carbon 
grains from Gail & Sedlmayr (1987). 

Due to the extreme temperature-dependence of the nu- 
cleation and dust evaporation process, the noise level in 
the temperature determination should be as small as ±1 K 
(see results in Sect. l3.1T)l . The statistical error in our MC 
method is determined by the number of crossing events 
of photon packages through a cell under consideration. 
Therefore, in order to achieve the required accuracy, one 
has to shoot many photon packages and one cannot use 
too small computational cells. 

This condition practically prevents a solution of the ra- 
diative transfer problem on the partly very highly resolved 
AMR hydro-grid. Instead, a comparably coarse regular 
grid (equally spaced in log r and 6) is used for the radia- 
tive transfer calculations. A gas and dust mass conserv- 
ing scheme has been implemented, based on cell volume 
overlaps, to map the {p, Ls}(r) AMR hydro-data on this 
transfer grid. Next, the opacities are calculated and the 
Monte Carlo experiment is carried out, which results in 
the mean intensity structure J(r) = -T r 4 d (r). This data 
is finally interpolated back onto the hydro grid, using a 
2D piecewise linear interpolation scheme. The whole pro- 
cedure is treated as additional source term in the frame 
of the operator splitting method of the hydro solver. 

In order to achieve the required temperature accuracy 
of ±1 K, the computational cost for one radiative transfer 
call is quite remarkable: about 40 sec on one processor for 
a 200 -grid (ID) and about 120 sec on 64 processors for a 
(120 x 80) -grid (2D). We are therefore forced to call the 
MC radiative transfer routine less frequent than at every 
hydro timestep. 

Instead, we have developed an extrapolation scheme 
that allows us to call the radiative transfer routine only 
at about every 3 rd to 100 th hydro timestep. The basic idea 
for this extrapolation in time is to interpolate in optical 
depth space. A parallel algorithm to compute the spherical 
optical depths 

-Rout 2 

TL(r,6,t) = J K ext (r',e,t) (^-j dr' (17) 

r 

on the basis of the AMR data has been implemented, sim- 
ilar to Rijkhorst et al. (2005). i? Q ut is the outer boundary 



ol the model domain and H" the stellar radius ol the ini- 
tial model (see Sect. I2.6[l . Given the last radiative transfer 
results T ro f(r) = T ra( j(r, t le { ) at reference time t le f with 
reference optical depths r re f(r) = 7L(r,t re f), the actual 
radiation temperatures are calculated as 

T r 4 ad (r, t) = T r 4 ef (r) + a cm (r)(t - t Ie{ ) 

+ &cor(T-) (r L (r,t) -T rof (r)) , (18) 

where a cor and 6 cor are local fit coefficients which are 
updated after every radiative transfer call. b COI (r) = 

d" d (r) is the local dependence of Tf ad from tl , which 

expresses typical backwarming effects. These effects can 
occur as fast as hydrodynamical changes, for instance if 
an opaque dust cloud enters a computational cell. 6 cor is 
calculated by linear regression of the local {T r 4 d , tl} data 
points along the radial direction. 

a cor expresses the net change caused by all other ra- 
diative transfer effects, e.g. by the passing of a dust cloud 
sideways of a considered cell or by the changing shadow of 
a growing dust cloud between the star and the cell. These 
changes are usually caused by the more distant opacity 
configurations and hence occur on much longer timescales. 
After having fixed 6 cor , the second fit coefficient a cor is cal- 
culated from inversion of Eq. Ijl8|l . such that the extrap- 
olation from the old to the new radiative transfer results 
would have been perfect. 

The time interval AtMC, after which the next proper 
MC radiative transfer must be calculated and the fit coef- 
ficients a cor and 6 cor are renewed, is controlled by the con- 
dition that the difference between the extrapolated Trad- 
values and the newly calculated T rc f MC results (measured 
by the 4-norm) must be smaller than 2 x AmcT, where 
AmcT ~ 1 K is the Monte Carlo noise in the temperature 
determination. 

We note that Eq. (|18J) allows for an exact treatment of 
the following two limiting cases 

1 4 3 4 

^rcf — q ^ cff ' acor = ^' cor = 4 "^ eff ' Trcf = ^ 
=>■ spherical diffusion approximation 
(Eddington approximation) 

T r 4 cf = W(r) 2* a cor = 0, b cm = \ T 4 (J^J , r rcf = 

=> spherical two-stream approximation 
(Lucy's approximation,) 

where W(r) = | [l — (l— 73-) is the radial dilution fac- 
tor. By fixing T ro f, a cor and b cor in this way, or by keep- 
ing them updated according to the actual values of the 
stellar radius i?*(t) and stellar temperature T+(t), respec- 
tively, one can run quick test models without involving the 
expensive MC radiative transfer routine. These simplifica- 
tions, however, are only possible for spherically symmetric 
ID models. 2D models require the proper inclusion of the 
MC radiative transfer. 

Equation l|18|) is hence a reasonable approximation 
with the ability to self-adapt to certain limiting cases if 



the ML results indicate so. In general, it will only be valid 
in a local sense for a limited time span after a last proper 
solution of the radiative transfer problem. The extrapo- 
lated radiative transfer results are always time-resolved 
because they depend on the actual optical depths. 

Regarding the MC radiative transfer routine itself, 
scattering and thermal re-emission are numerically treated 
in a different way (Niccolini, Woitke & Lopez 2003). 
However, re-emission and isotropic scattering is physically 
indistinguishable in the grey case, so we can treat the 
albedo 7 as free numerical parameter. Best performance 
is achieved by choosing 7 = 0.9. 

Another troublesome point is the Td-dependence of 
K^dust which makes necessary an iterative procedure be- 
tween opacity and radiative transfer calculations. We solve 
this problem by calling the MC routine only once (if it is 
foreseen for this timestep) and then using Eq. in the 
subsequent iterations. After typically 1-5 iterations, this 
procedure converges (opacity calculation - optical depths 
calculation - T ra d recalculation) . Only if it does not con- 
verge (which happens very seldom), more MC radiative 
transfer calls within this iteration loop are required. 

2.5. Boundary conditions 

Finding suitable boundary conditions for the multi- 
dimensional stellar wind problem in Eulerian coordinates 
(fixed volume) is not trivial. The main difficulty arises 
from the mass loss through the outer boundary which, 
on average, should be compensated for by a mass in- 
flow at the inner boundary. The FLASH-code uses 4 layers 
of guard cells to describe boundary conditions. On these 
guard cells, all variables (p, v, e, Lj) must be prescribed 
at every computational time step. 

2.5.1. Inner boundary condition 

The basic idea of our permeable inner boundary is to 
fix the gas pressure at the inner boundary. If matter is 
removed from the regions close to the inner boundary, 
e. g. blown away by radiation pressure, a pressure gradi- 
ent across the inner boundary develops which causes the 
desired mass inflow. 

In order to quantify this idea, we visit Eqs. (JIJ and J5J 
in the spherically symmetric stationary case, where they 
can be combined into the so-called wind equation 

dv r ( ci\ 2ci dcL 

V r ^— 1 5" h -7T- = arad ~ a grav . (19) 

or \ v£ J r or 

°t — \ppfp~ is the isothermal sound speed. Assuming 
v r <C ct (subsonic limit), Eq. (|19(l becomes 

1 dv r 2 | 1 gc| _ a rad - a grav 

v r dr r ci^ dr c\ 

By means of pv r r 2 = const (the stationary case of Eq.Q 
one can show that the l.h.s. of Eq. (|20J) equals In (pc|^), 



resulting m the integral 



z.b.5. Angular boundary conditions 



In (pc%) = ln(pc£) 



»rad 



■ dr 



(21) 



which shows that stationary flows in the subsonic limit 
in fact obey the hydrostatic condition. We use Eq. (|21|l 
to describe our inner boundary condition. We numerically 
evaluate the integral in fEa . I21j) from an arbitrary point r 
on the guard cells to the first valid cell within the model 
volume ro by assuming that the gas temperature T g can 
be linearly extrapolated backward from the first two valid 
cells in the radial direction. 

One way to fix the gas pressure at the inner bound- 
ary would be (/oc^)| r = p(ro,t — 0). However, this simple 
idea leads to the problem that a sudden increase of the 
temperature level close to the inner boundary (e. g. due 
to a dust shell formation event) leads to a likewise sud- 
den increase of the gas pressure level, which can result 
in a positive pressure gradient across the inner boundary, 
causing a vivid infall (negative velocities) through the in- 
ner boundary. According to our experience, it is better to 
fix the gas density at the inner boundary by 



{pel) 



p(r o ,t = 0) c^(r ,i) 



(22) 



which leads to a more regular increase of the pressure level 
around the inner boundary (also on the guard cells) in case 
of a sudden irradiation from the outside. A positive tem- 
perature offset lowers | 3 g" p |, i- e. it causes a modest out- 
ward acceleration of the gas in case of a sudden irradiation 
from the outside, which seems more natural regarding a 
star to be characterised by an equilibrium between gravity 
and pressure forces. 

Once having solved Eq. I|21|) for p, the radial velocity 
on the guard cells v r is calculated from pv r r 2 = const, the 
internal energy e is calculated from p and T g via the equa- 
tion of state, and the guard cells at the inner boundary 
are assumed to be dust- free £j(r<7- ) = 0. 

Since we do not consider stellar pulsation in this pa- 
per, a constant stellar luminosity L+ is assumed as inner 
boundary condition of the radiative transfer problem. 



2.5.2. Outer boundary condition 

For the outer boundary condition, a standard outflow con- 
dition (zero-gradient, i. e. a direct copy from the last valid 
cell designated by the subscript 0) is applied to the vari- 
ables v, T g and Lj. The gas pressure is extrapolated in 
such a way that the gradient compensates gravity and 
radiative acceleration (force- free) : 



In p 



= In po + -5- 1 GM, 
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(23) 



Equation (|2~3")l can be deduced from Eq. l(2"T|) under the as- 
sumption that M r , ct and K oxt are constant on the guard 
cells. Mass density p and internal energy e are finally cal- 
culated from p and T g , using the equation of state. 



In case of a 2D model, the guard cells at 8 < and 8 > ir are 
filled by a standard reflecting boundary condition which 
assures ve (r, 8 — 0, t) = ve(r,8 = TT,t) = 0. 

2.6. Initial condition 

The dynamical simulations are started from a dust-free 
static stellar atmosphere as initial condition. In the static 
(and hence spherically symmetric) case, the Poisson equa- 
tion, the equation of motion and the definition of the opti- 
cal depth form a system of ordinary differential equations 



dM r 
dr 
dp 
dr 
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(24) 
(25) 

(26) 



which is completed by the equation of state p — p(j>, T g ), 
the gas opacity K oxt = ^gas^Pi^g) an d an analytic for- 
mula of the temperature as function of optical depth 
T g = T g (r, tl). This formula represents the solution of 
the radiative transfer problem in radiative equilibrium. 

Starting with an initial guess of this dependency, 
Eqs. (|2"4*jl to (123) can be solved by a standard ODE solver 
with the outer boundary condition 



tl(Ro 







and two conditions at the stellar radius 



(27) 



(28) 
(29) 



where the initial stellar radius R® is given by the Stefan- 
Boltzmann law L+ = 47ri?° vT^. 

Having obtained this solution, the MC radiative trans- 
fer routine is called for given gas density distribution p(r) 
which provides new {r, T s ,tl} data points, from which 
an updated version of the temperature formula (piecewise 
linear interpolation) is calculated. The procedure needs 
about 3-10 iterations to converge. 

The physical parameters of the complete static stellar 
atmosphere problem are hence the stellar luminosity L+, 
the effective temperature T e s and the stellar mass M*. 
The static stratification is calculated during the initialisa- 
tion phase of the FLASH-hydrocode. The solution is then 
mapped onto the initial FLASH-grid, again using a piece- 
wise linear interpolation scheme in r. 

2. 7. Timestep control 

The computational timestep of the dust radiation hydro- 
dynamics code is limited by two criteria. First, 0.8xthe 
Courant-Friedrich-Levy timestep, given by the minimum 
of the sound wave travel time through any computational 



cell, and second, a temperature hmiter. Ihe latter limits 
the maximum change of T ra( j during one computational 
timestep to some given value. We choose this value to be 

5 K, which seems necessary for the simulation of dust for- 
mation and in fact stabilises the results. 

2.8. Model domain and grid resolution 

As radial model domain, the interval r/R® 6 [0.9, 10] is 
considered. In case of a 2D-model, the full angular range 

6 G [0, 7r] is considered. The inner boundary is chosen to 
be located at tl « 5, where the velocity variations are 
truly subsonic and the radiative transfer is diffusive. The 
outer boundary is located far outside the dust formation 
and wind acceleration zone. We have checked that a larger 
choice than 10 R® has only minor influence on the result- 
ing wind quantities like e. g. the mass loss rate. 

We start with a basic grid resolution (refinement 
level 1) of [rx6] = 96 x 64 cells. Based on second derivative 
criteria, the grid will be refined during the simulation to 
maximum level 5. Each refinement step increases the lo- 
cal resolution block-wise by a factor of two in each spatial 
dimension, such that the maximum resolution (if the grid 
needs to be fully refined) is 1536 x 1024. 

2.9. Computational limitations 

The implementation of the physical and chemical pro- 
cesses (source terms) in the frame of the FLASH-solver is 
independent of spatial dimension and applied coordinate 
system, such that now 1D/2D/3D Flash simulations with 
Cartesian or curvilinear coordinates are principally pos- 
sible. We find it, however, most convenient, if not even 
mandatory, to formulate the boundary conditions with 
spherical coordinates. Furthermore, the Monte Carlo ra- 
diative transfer code does currently only allow for ID and 
2D models with spherical coordinates. 

A complete high-resolution ID model (60 years of sim- 
ulation) requires about 5 x 10 5 CPU sec, a complete 2D 
high-resolution model about 5 x 10 7 CPU sec on Aster, 
which is an SGI Altix 3700 system consisting of 416 pro- 
cessors (Intel Itanium 2, 1.3 GHz). We have typically used 
8 processors for a ID model (16 hours user time) and 
64 processors for a 2D model (9 days). Hence, even with 
the present power of modern parallel super-computers, 3D 
simulations with Monte Carlo radiative transfer are simply 
too expensive for the required accuracy (see Sect. l3.l"T|l . 
Thus, for practical purposes, only ID (spherically symmet- 
ric) and 2D (axisymmetric) simulations can be run with 
the present code. 

3. Results 

The description of our results is separated into two parts. 
First, we compare our results for spherically symmetric 
(ID) dust-driven winds with other published works in 
Sect. 13.11 intending to validate our new numerical imple- 



lable 1. Kesults of ID spherically symmetric models 
for C/O increase (see text). Constant parameters: M* = 
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breeze: Subsonic solution with finite outer pressure for 
r—>oo (unphysical). In the numerical model, such 
solutions are featured by a standing shock wave in 
the outer regions. 

P = ... Quasi-periodic solution. All global quantities like 
M(t) or Md(t) show cyclic variations (see Fig. 0. 
The Fourier spectra are featured by an outstanding 
peak at a certain period P [days] . Dust is produced 
at every cycle, but small cycle-to-cycle variations 
may exist. 

irregular: Oscillatory solution with irregular dust shell 
production. No clear period. Sometimes period 
switches. 

Stationary solution with small oscillatory pertur- 
bations (very close to the bifurcation point). 

mentation. In Sect. l3~21 the new axisymmetric (2D) mod- 
els will be presented and discussed. 

3.1. ID results: code validation 

Fleischer et al. (1995) and Hofner et al. (1995) have pub- 
lished spherically symmetric (ID) models for time- 
dependent dust-driven winds of carbon stars with the 
same physical equations and assumptions (e. g. no pul- 
sation) and parameters (e. g. opacities) as in our model, 
but using different numerical techniques. These results can 
hence be used to validate our new computational imple- 
mentation. We will briefly summarise the basic results ob- 
tained by all three models, before we can discuss the dif- 
ferences. 

After the start of a new simulation, the initially hy- 
drostatic and dust-free stratification of the stellar atmo- 
sphere and the circumstellar environment undergoes some 
dramatic changes with strong shock waves in the outer re- 




Fig. 1. Stationary ID wind solution for M* = 1M , 
L* = 1O 4 L , T cff = 2400 K and C/O = 1.7. The thin 
black lines show the velocity v (full), the mass density p 
(dashed), the dust and gas temperatures Td and T g (full 
- practically identical in this model), the degree of con- 
densation / con d (full) and the nuclcation rate per hydro- 
gen nucleus J^/n^H) (dashed). The underlying thick grey 
lines show the results for a comparison model, where the 
Monte Carlo radiative transfer code is not called anymore 
(see text). The difference between black and grey indi- 
cates the effect of the Monte Carlo radiative transfer noise 
(ATmc ~ 1 K) on the solution. 



gions. However, the model structure soon relaxes towards 
a new type of solution, which is found to be either a sta- 
tionary wind or an oscillating wind. 

3.1.1. Stationary winds 

The stationary solutions are featured by low mass loss 
rates M, low outflow velocities and low dust-to-gas 
ratios Pd/Pg and occur for comparably small C/O ratios 
(see Table0. Figure^shows an example for such a finally 
almost time-independent solution. The inner regions, from 
the bottom of the atmosphere to the start of the dust 



Fig. 2. Snap shot of an ID oscillating wind solution for 
the same stellar parameters as in Fig. ^ but with slightly 
increased C/O = 1.75. In this plot, grey full lines refer to 
l.h.s. quantities and black dashed lines to r.h.s. quantities. 

formation zone, are approximately in hydrostatic equilib- 
rium. The density falls off exponentially here and the ve- 
locity is approximately zero. At some distance from the 
star (here at about 1.5 i?*) the temperature becomes low 
enough for nucleation which triggers the process of dust 
formation. Radiation pressure on dust then accelerates the 
dust /gas mixture and drives it through the sonic point lo- 
cated at about 2 i?* in the depicted model. The wind re- 
gion beyond that point is featured by an almost constant 
outflow velocity and a p oc r" 2 law. Due to the increas- 
ingly rapid outflow and the decreasing gas densities, the 
process of dust formation freezes in. For further details, 
see e. g. Gail & Sedlmayr (1987). 

Note that the formation of dust remains incomplete. 
Only about 18% of the condensable carbon has actually 
condensed into solid particles in the depicted model (see 
/cond in Fig.^l. The influence of the dust formation on the 
radiative transfer results is small. The model describes an 
optically thin, low-velocity, dust-driven stellar wind. 

Due to the Monte Carlo noise in the temperature de- 
termination of our radiative transfer code (see Sect. 12.4(1 . 
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Fig. 3. Examples for "quasi-periodic" (l.h.s.: T cff = 2500 K and C/0\1.6) and "irregular" (r.h.s.: T cff = 2400 K and 
C/0\1.6) ID oscillatory wind solutions. Depicted is the total gas+dust mass M(t) (upper row) and the total dust 
mass Md(t) (lower row) contained in the model volume. The change of M(t) is mainly due to mass inflow and outflow 
through the inner boundary. The change of M&(€) is due to dust formation and mass loss through the outer boundary. 
In the quasi-periodic model (l.h.s.), a new shell of dust is formed at every second oscillation cycle (double-periodic 
solution). Fourier analysis shows two sharp maxima at 1060 and 530 days. The basic period P — 530 days given in 
Table prefers to the oscillation of the stellar atmosphere. The irregular model (r.h.s.) shows the production of dust 
shells at almost unpredictable instants of time. 



the relaxation towards a stationary wind is actually never 
complete. The slightly random radiative transfer results 
cause small temperature changes in the inner regions 
which create inward and outward travelling sound waves. 
The outward travelling waves amplify in the large density 
gradient close to the star. This amplification of waves and 
the strong temperature dependence of the nucleation pro- 
cess result in noticeable variations of the fluid variables in 
the outer regions of order 5%, although the temperature 
noise ATmc is as small as 0.1% (not visible). This result 
leads us to the conclusion that an accuracy of the radia- 
tive transfer results of ATmc ~ 1 K is actually required 
to model dust-driven winds, which makes the simulations 
very expensive. We have checked that a reduction of the 
Monte Carlo noise by a factor of two (by using four times 
more photon packages) reduces the amplitudes of the re- 
maining disturbances by about a factor of two. 

We can suppress the creation of new disturbances com- 
pletely by not calling the radiative transfer routine any- 
more, thereby fixing T re f and 6 cor for ever and putting 
a cor = (see Ea, ITH|) 2 . In that case, the stationary solution 
becomes completely smooth and freezes in for all times 
(see underlying grey model in Fig.^). This freezing for ar- 
bitrary long times is possible, because our inner boundary 
condition compensates the mass loss through the outer 
boundary, in contrast to the Fleischer et al. (1995) and 
Hofner et al. (1995) models, where the model volume runs 
empty for very long times. 

3.1.2. Oscillatory winds 

Beyond some critical value for the stellar luminosity 
(Fleischer et al. 1995) or the C/O ratio (Hofner et al. 
1995), the stationary wind solutions become unstable. 
This instability, called "exterior K-mechanism" , is caused 
by feedbacks of the dust formation process on the radia- 
tive transfer. After the formation of a new optically thick 



2 This is not the same as fixing the temperature profile: t- 
changes still cause T ra d-changes. 



dust shell, the temperatures inside of this shell increase 
quickly due to backwarming (see 2d(r) in Fig. 0), which 
switches off the nucleation process and may even cause 
dust evaporation in the inner regions. The newly formed 
dust shell is hence truncated at its inner edge and driven 
outward, which reduces again its optical depth because of 
radial dilution. Consequently, the temperatures in the in- 
ner regions start to slowly decrease again, and a new cycle 
of dust formation may begin. 

This mechanism, however, interferes with other hydro- 
dynamical processes in the dust formation zone. The sud- 
den acceleration of a newly formed dust shell creates two 
waves. The first wave, an outward travelling compression 
wave, steepens up into a shock wave, compresses the gas 
ahead of the outward moving shell and facilitates sub- 
sequent dust formation in the wake of this shock. The 
second wave, an inward travelling expansion wave, causes 
the dust formation zone to be re-filled with fresh matter 
from the star. The sudden increase of the temperatures 
in the stellar atmosphere creates a third wave, a pressure- 
driven expansion wave, which merges with the incoming 
expansion wave to create a new shock wave below the dust 
formation zone. These hydrodynamical disturbances affect 
the structure of the stellar atmosphere and may interfere 
with the stellar pulsation (disregarded in this paper). The 
characteristic timescale of the dust-induced K-mechanism 
(usually longer than typical pulsational periods) may or 
may not coincide with the timescale of these hydrodynam- 
ical processes. A complicated time-dependent behaviour 
results in this way, which can be periodic, multi-periodic 
or completely irregular (see Fig.[3J). 

As a common result from all three simulations (this 
work, Fleischer et al. 1995, Hofner et al. 1995), the oscil- 
latory solutions are featured by higher mass loss rates 
(M), higher outflow velocities (voo) and higher dust-to- 
gas ratios (pd/p g ) as compared to the stationary solutions 
(measured by temporal mean values at the outer bound- 
ary). The solutions describe optically thick winds of some- 
times massively obscured central stars, e. g. infrared car- 
bon stars. 



Ihe effect of the Monte Carlo noise on the oscilla- 
tory wind solutions is more difficult to assess, because 
we cannot suppress it completely. The ongoing creation 
of small temperature disturbances by the MC noise in 
the stellar atmosphere prevents a complete relaxation to 
a truly periodic solution. Therefore, our solutions are al- 
ways featured by small deviations from periodicity or are 
even irregular. This behaviour, however, was also found by 
Fleischer et al. (1995) and Hofner et al. (1995), who used 
deterministic radiative transfer methods. If we use the de- 
terministic Lucy approximation for the radiative transfer 
(see Sect. we can produce solutions for some param- 
eter combinations which are sometimes closer to truly pe- 
riodic, but the general behaviour of the solutions is very 
similar. Therefore, we conclude that the Monte Carlo noise 
is not a critical issue concerning the resulting type of wind 
solution. 

3.1.3. Comparison to other ID models 

The general behaviour of our new Flash simulations for 
dust-driven C-star winds in the special case of spheri- 
cal symmetry resembles well the behaviour reported by 
Fleischer et al. (1995) and Hofner et al. (1995), progress- 
ing from breezes — * stationary winds — > oscillatory so- 
lutions with increasing C/O. In a comparative study, 
Hofner et al. (1996) have argued for good general agree- 
ment between the two codes despite the different numer- 
ical techniques. Since Fleischer et al. (1995) have consid- 
ered throughout higher stellar luminosities as in this pa- 
per, we will concentrate on the comparison to the Hofner 
et al. models in the following. 

Figure 0] shows the differences of the calculated mean 
wind quantities to Hofner et al. (1995, see their Table 3). 
In summary, the Hofner et al. winds are a little stronger 
with slightly larger M and . One can also describe these 
differences in the following way. The "jump" from sta- 
tionary — » oscillatory solutions in our models occurs at 
a slightly larger critical C/O value as compared to the 
Hofner et al. models. The obtained periods of Hofner et 
al. show the same trend (increasing with increasing T c ff) 
but are throughout shorter by a factor of about 1.5. 

The slightly lower efficiency of the driving by dust in 
our work can already be explained by a tiny detail in 
the treatment of equilibrium chemistry. In contrast to 
Hofner et al. (1995), we have included the C3 molecule 
which reaches high concentrations around « 2000 K and 
thereby reduces the number density of C atoms and C2H2 
molecules. This lowers the efficiency of the nucleation and 
growth of the carbon dust particles. According to our 
model, just ignoring the C3 molecule results in higher M 
and Voc by a factor of about two in the stationary regime, 
which would put our results even above those of Hofner 
et al. in Fig. ||l 

However, various other numerical details could also be 
the reason for the shown discrepancies. With regard to the 
very complex physical system under investigation, small 
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Fig. 4. Comparison of ID results to Hofner et al. (1995) 
for M* = 1 Mq , L± = 1O 4 L and T cff = 2500K. 



Table 2. Results of ID spherically symmetric models 
for C/O decrease (see text). Constant parameters: M* = 
1M© and L„ = 10 4 L Q . 
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osc. breeze: Low- velocity solution (sometimes subsonic Voo) 
with small oscillations and persistent dust inside 
the nucleation regime (see Fig. |HJ . 

^ : Double-periodic: formation of a dust shell at ev- 

ery second oscillation cycle of the outer stellar 
atmosphere. 



numerical deviations can have big effects. In view of the 
different radiative transfer methods, the different imple- 
mentation of the hydrodynamical equations with differ- 
ent numerical schemes (implicit/explicit) and, last but not 
least, the different inner boundary condition, the degree 
of agreement is actually quite remarkable. 

3.1.4. The folding bifurcation 

The results of our ID simulations, as summarised in 
Tabled] have been obtained by a stepwise increase of C/O, 



O 



breeze j 



station^ 



C/O 

Fig. 5. Folding Bifurcation in modelling C-star winds. £7 



giving the model sufficient time in between to relax to the 
respective equilibrated solution (about 50yrs). 

However, if we stepwise decrease C/O again, we face 
a surprise (see Table EJl ■ The model does not fall back 
toward the stationary solution at the critical C/O value 
for increase as observed before. Instead, it stays oscilla- 
tory until it finally falls back to a "breeze" at another, 
much lower critical C/O value for decrease. Apparently, if 
once the transition to an oscillatory mode is made and 
the outer parts of the stellar atmosphere are swinging 
back and forth considerably, a more efficient mechanism of 
dust production can be sustained even with less fuel. This 
non- uniqueness for intermediate C/O values (hysteresis 
behaviour / folding bifurcation) is sketched in Fig. and 
has not been reported so far. 

3.1.5. Oscillating low-velocity solutions 

The reason for the folding bifurcation sketched in Fig. [3] 
lies in an accumulation of small amounts of dust particles 
in the "metastable" region close to the star which is too 
hot for nucleation, but too cold for dust evaporation. In 
a stationary wind, this region remains completely dust- 
free because of the throughout positive radial velocities. 
However, in case of oscillations, some dust particles (which 
have been created farther out) periodically re-enter this 
region, where they temporarily stay and grow up to very 

1/3 

big particles with radii as large as (a) = (J^) L\/Lq » 
50 /im. The total surface of these few big grains remains 
so small, however, that an average molecule needs about 
Tcoii = (t , th(367r) 1 / 3 / oL 2 ) ~ 100 yrs to collide with them. 
Hence, the degree of condensation and the dust opacity 
are practically frozen in and remain too small to cause a 
net outward acceleration |a ra( j| Ss |a grav |. However, a rac j is 
just large enough to change the pressure gradient in this 
quasi-static region, which leads to higher densities in the 
outer parts. 

This density levitation facilitates dust formation in the 
outer parts and leads to a persistent, oscillatory stellar 
outflow with low (sometimes subsonic) velocities and a 
smooth dust distribution, see Fig. El Winters et al. (2002) 
described similar low- velocity, almost quasi-static ("B- 




21 ° 



Fig. 6. Subsonic ID solution with small oscillations for 
Ah — 1 Mq , L± = 1O 4 L , T cff = 2600 K obtained by de- 
creasing C/O \ 1.6 after about 1800 years of simulation. 
Note the small but non-zero degree of condensation inside 
the nucleation peak. For clarity, we have "frozen" this so- 
lution lately by not calling the MC radiative transfer rou- 
tine anymore (see Fig. ^ and text). 



type") ID wind solutions in a certain range of stellar pa- 
rameters under the influence of a strong stellar pulsation 
as inner boundary condition. Here, we find a same type of 
wind solution even without stellar pulsation. 

The question remains open whether or not such low- 
velocity outflows exist in reality. Jura et al. (2002) noted 
that the close-by semi-regular variable L 2 Pup (spectral 
type M5 IIIe-M6 Hie) has in fact a very low- velocity out- 
flow (voo » 1.9km/s). Based on continuum observations, 
Jura et al. (2002) estimated the mass-loss rate of L 2 Pup 
to be M w 3.5 x 10 -7 M Q /yr, whereas based on molecu- 
lar line observations, Winters et al. (2002) estimated Ms=s 
2 x lO _9 M0/yr. Despite this large scatter, these wind 
properties seem to match the characteristics of the low- 
velocity outflows as found by Winters et al. (2002) and in 
this work. However, L 2 Pup is an oxygen-rich star with 
lower luminosity, higher stellar mass and higher effective 
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Fig. 7. Rayleigh- Taylor instabilities in an expanding dust 
shell shortly after the start (t = 3.2 yrs) of a new axisym- 
metric (2D) simulation. Stellar parameters M* = 1M@, 
L+ = 1O 4 L , T cS = 2500K and C/0 = 1.9. The additional 
contour line for T ra( j = 2500 K indicates the size of the star. 

temperature as assumed in these models, so this link be- 
tween observations and theoretical models is thin. 

3.2. Results of the 2D models 

The major aim of this work is to advance from ID to 
2D dust-driven wind models. This enables us, for the first 
time, to explore the effects of deviations from spherical 
symmetry on dust-driven winds in a consistent way 3 . Can 
these complicate physico-chemical systems, which have 
been proven to be unstable already in ID models, keep 
their spherical symmetry, or — do instabilities lead to 
spontaneous symmetry breaking, structure formation and 
more complex flow patterns, for example turbulent mix- 
ing? The answer of the 2D simulations is quite unambigu- 
ous. 

3.2.1. The Rayleigh-Taylor type instability 

Already shortly after the start of a new simulation, the 
expanding dust shells experience Rayleigh-Taylor type in- 
stabilities at their outer edges (see Fig. 01, which break the 
spherical symmetry spontaneously. This instability will be 
further analysed and discussed in a forthcoming paper 
(Woitke, Farzinnia & Weterings 2006), but their nature 
is easy to understand. The dust containing gas, which has 
a larger opacity K oxt and hence a higher potential energy 
— GM t (l - a)/r with a = K cxt L*/(47rcGM*) (see Eq. EJ), 
is accelerated into a dust-poor gas with lower opacity 

3 Freytag & Homer (2003) have passively integrated the dust 
moment equations in first 3D simulations, i. e. without the im- 
portant feedbacks of the dust on the acceleration and on the 
radiative transfer. 



and hence lower potcntiaf energy, thus, it is cnergcticalfy 
favourable to exchange matter across the outer interface 
of a dust shell, which makes the fluid unstable. 

According to the 2D simulations, a spontaneous mix- 
ing occurs in most cases at the outer edges of dust shells, 
where dust-rich fingers penetrate outward and more trans- 
parent dust-poor fingers penetrate inward with respect to 
the outward moving shell. The occurrence of Rayleigh- 
Taylor type instabilities is a persistent feature also in the 
more relaxed 2D models. 

3.2.2. Description of 2D solutions 

In the following, we will describe the qualitative behaviour 
of two relaxed 2D models with emphasis on the main new 
features caused by dropping the assumption of spherical 
symmetry 4 . 

The first high-resolution 2D model has been calculated 
for parameters M* = 1 M Q , L± = 10 4 L Q , T cff = 2500 K 
and C/O = 1.9. This set of parameters results in an op- 
tically thin stationary wind with only little ^-variations. 
The radial structures (cuts along constant 9) resembles 
well the respective ID model structure (see Sect. 13.1. Ti l. 

The second high-resolution 2D model has been cal- 
culated for the same stellar parameters, but now with 
C/O = 2.2, a parameter combination which resulted in an 
optically thick wind with irregular oscillating behaviour 
in ID (see Table 0). Figure [S] shows some details of this 
2D model after 46.5 yrs of simulation, after which the re- 
spective ID model is sufficiently relaxed. 

The most striking new feature is that the dust "shells" 
are not complete, but do only occupy a certain frac- 
tion of the total solid angle. These dust "arcs", some- 
times also smaller "caps" originate from the dust forma- 
tion zone which is mainly radially oscillating according to 
the exterior K-mechanism known from the ID models (see 
Sect. I3.L2*|) . However, these oscillation often become out- 
of-phase concerning north-polar, equatorial and south- 
polar regions, etc., in the final relaxed state. Consequently, 
dust forms from time to time in restricted regions above 
the stellar surface. Only in those cases where the complete 
dust formation zone around the star performs a concerted 
radial oscillation, a complete radial dust shell is produced. 
In all other cases, the production of limited dust arcs is 
more typical. 

In fact, the different polar and equatorial regions seem 
to perform more or less their own independent oscilla- 
tions, since they are only weakly coupled to each other 
via tangential velocities and radiative transfer feedbacks. 
Considering the often slightly chaotic behaviour of the 
dust formation zone in the ID models (weak chaos), 
it's clear that small disturbances (e. g. introduced by the 
Monte Carlo noise or by the Rayleigh-Taylor instabil- 

4 Unfortunately, the 2D models are so expensive (see 
Sect. |2~§! that only two complete high-resolution 2D simula- 
tions are available up to now, supplemented by a number of 
additional low-resolution simulations. 
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Fig. 8. Snap-shots after 46.5 years of simulation from an axisymmetric (2D) model with parameters M* = 1M Q , 
L+ = 10 4 L Q , T c g = 2500K and C/O = 2.2, showing the following quantities as grey-scale contour plots: first row: 
mass density log p and radiation temperature T ra d, second row: degree of condensation / con d and nucleation rate 
per hydrogen nucleus log(J+/ri(H)), last row: radial velocity v r and tangential velocity vg (white indicates a motion 
towards the north pole and black towards the south pole). The additional black contour line for T ra d = 2500 K indicates 
the size of the star. Central blank regions are not included in the model. 



rty) may cause completely different trajectories alter sul- 
ficiently long times. 

There is, however, one mechanism that tends to equili- 
brate the phases of the oscillations and to reinstall spher- 
ical symmetry, namely the radiative backwarming, which 
interrupts the dust production inside a forming shell (see 
Sect. l3~l~2")l . However, this backwarming is not as simple 
as in ID geometry. It has only a limited reach and cannot 
trigger the dust formation at the opposite side of the star. 
We believe that this is the reason why the typical angular 
extent of the dust arcs in the model is about one forth to 
one half of the total solid angle: An arc of such an angular 
extension produces locally almost as much backwarming 
as a complete radial shell. 

Figure El shows further details of the calculated tem- 
perature structures. The radiative transfer through the 
clumpy dust arcs causes complex temperature patterns 
with backwarming, sidewarming and shadow formation. 
The consequence is that the next dust formation episode 
will occur preferentially in those regions close to the star 
which are now not so much affected by the radiative back- 
warming. 

Once a dust arc has formed, it is accelerated outward 
by radiation pressure. This outward motion causes not 
only a radial displacement and dilution, but also a tan- 
gential expansion because of the decreasing ambient pres- 
sure. Supersonic tangential velocities result in this way 
(see lower r.h.s. part of Fig. [HJ which increase the angu- 
lar size of the arcs. At the side edges of the dust arcs, 
the matter stays behind and continues to form more dust 
which leads to an apparent inward bending and extension 
of the arcs. Concerning the smaller dust caps, the expan- 
sion motion resembles the behaviour of mushroom clouds: 
tangential expansion at the top and counter-rotating vor- 
tices in the wake with Kelvin-Hclmholtz instabilities. 

The occurrence of negative (infall) velocities under- 
neath and within the dust formation zone is a typical fea- 
tures in all oscillating wind models, not only in ID but 
also in 2D. However, in the 2D case, dust-poor gas may 
fall back toward the star at the same time as dust-rich 
gas is accelerated outward. In fact, these infall motions 
occur typically just between the outward moving arcs and 
clouds. Since the infalling gas is compressed by the am- 
bient pressure, the infalling dust-poor regions are usually 
much smaller than the outward moving dust-rich struc- 
tures. 

The sudden radiative acceleration of the newly formed 
dusty regions leaves behind temporary gaps in the dust 
formation zone which are not only refilled by dust-free 
matter streaming in from the stellar atmosphere, but also 
by infalling matter falling back from the outside. These 
two streams collide in the dust formation zone, which cre- 
ates turbulence. This turbulence provides additional rea- 
sons for the subsequent formation of inhomogeneous dust 
distributions (Helling et al. 2004). 

Furthermore, the infalling matter may already contain 
some small amounts of dust particles which have been cre- 
ated by nucleation farther out. In this way, dust particles 




2.0x10" 4.0x10" 6.0x10" 8.0x10" 1.0X10" 1.2x10 u 



Fig. 9. Details of the radiation temperature structure 
in the 2D model shown in Fig. [5] (same timestep) with 
overplotted white contour lines for T ra d = 1000 K, 1200 K 
and 1500 K. The additional black contour line for K cxt = 
5 x 10 cm -1 encircles the optically thick dust struc- 
tures. Note the backwarming of the dust arc at the bottom 
and the formation of shadows behind the opaque regions. 

Table 3. Comparison of calculated wind quantities be- 
tween ID and 2D models for M* =lM 0) L t = 10 4 L©, 
T eS = 2500 K and C/O = 2.2. 



geometry 


(M) [Mo/yr] 


(voo)[km/s] 




ID (spheric, sym.) 


7.9(-6) 


31.7 


3.2(-3) 


2D (axisymmetric) 


5.6(-6) 


29.2 


3.2(-3) 



reach again the much denser dust formation zone, where 
they will serve as additional nucleation seeds. This pro- 
cess results to be the main cause for new dust formation 
events. Considering the history of sample grains, we find 
typically 

nucleation — > infall — > growth — > acceleration. 

Last but not least, the more distant dust shells and 
arcs exhibit a variety of small-scale sub-structures. These 
small-scale structures are the result of the action of vari- 
ous flow instabilities, not only the Rayleigh- Taylor insta- 
bility as briefly explained in Sect. 13.2.11 These instabili- 
ties apparently need some time to amplify perturbations 
and so to further shape and fragment the outward mov- 
ing dust shells and arcs. The result are cauliflower-like 
shapes, which could also be described as numerous small- 
scale clouds within the expanding dust shells and arcs. 

3.2.3. Comparing ID with 2D results 

Considering the calculated mean wind quantities (e. g. the 
mass loss rate), the results of the ID and 2D models are 



remarkably similar (sec lablcyj, despite all the detailed 
effects described in Sect. 13.2.21 The slightly lower mass 
loss rate in the 2D model (30%) could indicate that dust 
arcs are slighly less effective in lifting the dust-poor gas 
in between out of the gravitational potential of the star 
than complete dust shells. However, this deviation could 
also be a simple consequence of the insufficient time given 
to the 2D model to relax. 

The deviations along one specific line of sight, however, 
may be quite different. Consequently, the spectral appear- 
ance of the star can be strongly affected by the deviations 
from spherical symmetry presented in this paper. This will 
concern lightcurves, images and visibilities, in particular 
at short (optical and near infrared) wavelengths, where 
the flux consists mainly of dust attenuated stellar pho- 
tons. These effects need to be further investigated in a 
forthcoming paper. 

4. Summary and discussion 

New axisymmetric (2D) hydrodynamical models for 
carbon-rich AGB star winds have been developed which 
include coupled grey continuum Monte Carlo radiative 
transfer, equilibrium chemistry, time-dependent dust for- 
mation and radiation pressure on dust. According to our 
knowledge, this is actually the first successful application 
of the Monte Carlo radiative transfer technique in the 
frame of multi-dimensional hydrodynamical simulations. 

The spherical symmetric (ID) version of the code has 
been checked against other published works (Fleischer 
et al. 1995, Hofner et al. 1995) showing good agreement 
in qualitative wind behaviour and quantitative results. 
Considering the most simple case without stellar pulsa- 
tion (hydrostatic inner boundary condition), the ID sim- 
ulations reveal 3 basic types of dust-driven stellar winds: 
stationary winds with low mass-loss rates and low but su- 
personic outflow velocities, oscillating winds with shock 
waves and dust shells, high mass-loss rates and high out- 
flow velocities, and oscillating low-velocity outflows with 
low mass-loss rates and sometimes even subsonic outflow 
velocities. 

All three types of dust-driven winds may occur for the 
same stellar parameters M*, L+ and T c g if the carbon-to- 
oxygen ratio C/O is considered as free parameter. The de- 
pendence of wind type (and resultant wind properties) on 
C/O is found to be non-unique, i. e. dependent on history. 
By increasing C/O, we find a sudden jump from station- 
ary to oscillating wind solutions, whereas for decreasing 
C/O, we find a sudden jump from oscillating winds to 
oscillating low- velocity outflows. This folding bifurcation 
(hysteresis behaviour) has not been reported so far. The 
third type of wind solution (the oscillating low-velocity 
outflows) seem to resemble the B-type wind solutions pro- 
posed by Winters et al. (2002), although we have not con- 
sidered any stellar pulsation in this paper. 

The possibility of a star with given parameters M*, 
, T c ff and C/O to produce two different types of winds 
with high and low mass-loss rates might open up the pos- 



sibility of spontaneous or induced switches between these 
two mass-loss modes, which could be related to the ob- 
servation of distant radial dust shells (see e. g. Mauron 
& Huggins 1999). A mechanism to obtain such switches 
has been reported by Simis et al. (2001). 

The results of the new axisymmetric (2D) dust-driven 
wind models reveals an even more complex picture of the 
dust and wind formation around carbon-rich AGB stars. 
Although the calculated mean wind properties like (M) 
and (woo) of the new 2D wind models are rather similar 
to those of the ID models, the details exhibit several new 
effects. 

Despite the spherically symmetric initial and bound- 
ary conditions, spontaneous symmetry breaking occurs in 
these winds due to action of various instabilities. Dust 
is found to form from time to time in restricted regions 
above the stellar surface. These spatially restricted dust 
formation events lead to the production of incomplete dust 
shells, dust arcs or even smaller dust caps. Only in those 
cases where the complete dust formation zone around the 
star performs a concerted radial oscillation, a complete 
dust shell is produced. 

The opaque but geometrically thin dust structures are 
accelerated outward by radiation pressure, expanding ra- 
dially and tangentially like mushroom clouds, which leads 
to an accumulation of overlapping arcs in the more dis- 
tant regions. The optically thin, dust-poor matter in be- 
tween tends to fall back toward the star, not only temporal 
(between the active dust formation phases) but also spa- 
tial (between the outward moving dusty structures). This 
matter mixes with fresh dust-free matter from the stellar 
atmosphere in a turbulent way in the dust formation zone, 
which again leads to a subsequent production of irregular 
dust structures. 

Further away from the star, flow instabilities (e. g. 
Rayleigh- Taylor) have time to modify and fragment the 
outward moving dust arcs and shells, producing numer- 
ous small-scale clumps in the outward moving dust shells 
and arcs. 

From the few models calculated so far, we conclude 
that whenever the ID model relaxes to an irregular oscil- 
latory wind solution, the 2D model is likely to show strong 
deviations from spherical symmetry. This type of solution 
is typical for optically thick winds of embedded IR car- 
bon stars. In contrast, if the ID model relaxes toward a 
stationary wind solution (which is the typical result for 
optical carbon stars with optically thin winds), the 2D 
model shows a similar behaviour. 

The formation of incomplete shells, arcs or smaller 
caps due to an irregular dust production of the central 
star may be important to understand recent interferomet- 
ric observations of embedded carbon stars and red super- 
giants (e. g. Monnier et al. 2004), which show obvious de- 
viations from spherical symmetry. Woitke & Quirrenbach 
(2005) have calculated IR images and visibilities from the 



calculated ZD models discussed m this paper J which can 
be useful for the interpretation of such observations. 
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